Consequences of local gauge symmetry in empirical tight-binding theory 
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A method for incorporating electromagnetic fields into empirical tight-binding theory is derived 
from the principle of local gauge symmetry. Gauge invariance is shown to be incompatible with 
empirical tight-binding theory unless a representation exists in which the coordinate operator is 
diagonal. The present approach takes this basis as fundamental and uses group theory to construct 
symmetrized linear combinations of discrete coordinate eigenkets. This produces orthogonal atomic- 
like "orbitals" that may be used as a tight-binding basis. The coordinate matrix in the latter basis 
includes intra-atomic matrix elements between different orbitals on the same atom. Lattice gauge 
theory is then used to define discrete electromagnetic fields and their interaction with electrons. 
Local gauge symmetry is shown to impose strong restrictions limiting the range of the Hamiltonian 
in the coordinate basis. The theory is applied to the semiconductors Ge and Si, for which it is shown 
that a basis of 15 orbitals per atom provides a satisfactory description of the valence bands and 
the lowest conduction bands. Calculations of the dielectric function demonstrate that this model 
yields an accurate joint density of states, but underestimates the oscillator strength by about 20% 
in comparison to a nonlocal empirical pseudopotential calculation. 
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I. INTRODUCTION 

Tight-binding theory was originally proposed as an 
ah initio technique for calculating the electronic propt. 
erties of crystalline solids from atomic wave functions.u 
However, first-principles calculations based on a linear 
combination of atomic orbitals (LCAO) are computa- 
tionally very demanding, and the tight-binding approach 
met with relatively little success until Slater and Koste* 
suggested that it be used as an interpolation scheme,El 
in which the Hamiltonian matrix elements are fitted 
to experimental data or to band structures computed 
by other methods. This made it possible to describe 
atomic-level physics in a basis of minimal size, leading 
to wide-ranginaYM|H|)lications in many areas of condensed- 
matter physics Jj'clErElLl'El With modern computer capabil- 
ities, first-principles electronic-structure calculations are 
now commonplaca, and ah initio tight-binding theories 
are flourishing Jj'QS Yet even today, the empirical theorya 
predominates (even for the fitting of first-principles cal- 
culations) because it is simple and physically intuitive. 

The formalism of Slater and Kosterla is incomplete, 
however, in that it contains no prescription for coupling 
the electronic systenL.to_external electromagnetic fields. 
In ah initio theoriesJSuu one can use minimal ccwpling 
(with suitable modifications for nonlocal potentialgj) and 
calculate directly the necessary matrix elements of the 
momentum or velocity operator. In the empirical theory, 
these matrix elements can simply be treated as extra fit- 
ting parameters ,£30113 determined by fitting the dielec- 
tric function (and thus oscillator strengths) to experi- 
mental or first-principles spectra. However, even with 
the full use of symmetry restrictions, the number of addi- 
tional parameters can be undesirably large; for example, 
Chang and AspnesM have proposed an sp^(P model for 
GaAs with 13 Hamiltonian parameters and 17 indepen- 
dent momentum parameters. 

It is therefore clearly desirable to find ways of reducing 



or eliminating these extra fitting parameters. One possi- 
bility is to define a kinematic momentum operator (equal 
to mass m times velocity) by 



(1.1) 



where x is the coordinate of the electron and H is the 
Hamiltonian. In a sense this merely trades one prob- 
lem for another, since the coordinate matrix elements are 
still unknown, and the number of these allowed by sym- 
metry is no less than the number of momentum matrix 
elements. However, it is physically reasonable to simplify 
the coordinate matrix by setting all nonlocal matrix ele- 
ments to zero: 



(1.2) 



Here \a, x^) is the ket vedjoc-for an orthogonalizcd atomic 
orbital (Lowdin orbitallj'Efl) of type a located at posi- 
tion Xi. The parameter Xq,o,/ is an intra-atomic matrix 
element coupling orbitals a and a' on the same atom. 
The simplest choice of all is to set 



0; 



in this model, there are no fitting 



(1.3) 



those found in the HamihonianffijSffifi'BBS 
A closely _p^plfl±cd_._approach is the Peierls 
substitution,ElE3E3'y'BEl in which the zero-field 
Hamiltonian matrix (a, Xi|_ff x^/) for a particle of 
charge e is replaced by 



{a,yLi\H\a' ,yie) exp 



he 



A • dx + e(j){yii)5u>8o 



(1.4) 



in the presence of a vector potential A and scalar poten- 
tial 0. If the path oLptegration in Eq. (1.4) is chosen 
to be a straight line,B'C§Ea then the linear term in the 
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Taylor series expansion of this equatio n is t he s ame as 
the A • p coupling obtained from Eqs. (l.l)-([l^).C^ 

The total elimination of extra fitting parameters makes 
this model an attractive one. However, by eliminating 
the intra-atomic matrix elements Xq,q,/, one obtains a 
tight-binding model that is not valid in the tight-binding 
limit of isolated atoms. Thus, although the model should 
provide a reasonable description of inter-atomic transi- 
tions between extended states, one has less confidence 
in its ability to describe localized states, which may 
be important at surfaces or interfaces. Many authors 
have therefore suggested augmenting the zero-parameter 
model by including a. small. number of intra-atomic ma- 
trix elementsHZi'mme'El It has been shown for 
porous Si that, although the intra-atomic matrix el- 
emeiits are small in magnitude (in a Bloch-function 
basiad), the interference between these terms and the 
inter-atomipjpatrix elements contributes 25% of the total 
absorption.ca Thus, it appears that a quantitative treat- 
ment of nanostructures may not be possible (in general) 
without the inclusion of intra-atomic mat rix elenjents 

The main difhculty with such modcl£jWBElEiy'El 
is that they are not gauge invariant .E3 As shown by the 
examples in Refs. |^ and lack of gauge invariance can 
lead to gross qualitative errors in the predicted values of 
physical quantities. Thus, there are significant problems 
with both approaches considered above. The models with 
Tiaa' — are gauge invariant, but they cannot describe 
intra-atomic transitions. The models with Xqq/ 7^ can 
describe intra-atomic transitions, but they are not gauge 
invariant. 

The purpose of this paper is to demonstrate a tech- 
nique for constructing tight-binding models that are 
gauge invariant and provide a full description of intra- 
atomic transitions. This is achieved by treating empir- 
ical tight-binding theory not as an approximation de- 
rived from the Schrodinger equation, but as a fundamen- 
tal quantum- mechanical system in its own right. This 
theory is required to satisfy all of the basic principles 
of quantum mechanics, the most important of which 
(in the pr esgn ipCpa|tcxtj) is the principle of local gauge 
symmetryJ3S'LllL3[!3o The essence of this principle is 
the concept that electromagnetism in quantum mechan- 
ics is the gauge-invariant manifestatkw-ef a nonintegrable 
(i.e., path-dependent) phase factor 

As will be shown in Sec. || below, the 
isting models with intra-atomic coupling 
are not gauge invaciant is that the coordinates x, y, and 
z do not commute.c3 Gauge invariance requires commut- 
ing coordinates, the existence of which implies the exis- 
tence of a basis of coordinate eigenkets. Since empiri- 
cal tight- binding theory deals with finite vector spaces, 
the coordinate basis is necessarily discrete. Hence, the 
most general gauge-invariant finite vector space is a 
set of discrete coordinate eigenkets. This basis may 
be transformed to a tight-binding basis by constructing 
"orbitals" from symmetrized combinations of coordinate 
eigenkets (using well-known techniques for symmetrizing 



plane wave: 



The concept of gauge symmetry on a discrete lat- 
tice is not new, having appeared many years ago as a 
technique for imnoaing .a momentum cutoff in quantuia 
chromodynamics.E3E3E3u3c3 Governale and UngarelliEa 
have recently suggested using lattice gauge techniques in 
empirical tight-binding theory. However, their proposal, 
like most applications of lattice gauge theory, is based on 
a simple cubic lattice. As shown below, the simple cubic 
lattice is unsuitable for practical tight-binding models be- 
cause it can only achieve sufficient accuracy with an un- 
reasonably large basis (i.e., a very small lattice constant). 
Thus, the development of efficient tight-binding models 
requires consideration of more general geometries. 

Christ, Friedberg, and Lee0'@'@ have developed a 
lattice gauge theory for random lattices, which (with 
some slight modifications) is sufficiently general for the 
present purposes. However, the complete formal machin- 
ery of quantum chromodynamics is somewhat cumber- 
some when one is dealing only with simple electromag- 
netism. Thus, for reasons of clarity, the author has cho- 
sen to present the theory terms of a simple but elegant 
approach used by DiracE3 After a preliminary discus- 
sion of topology (i.e., how an electron is permitted to 
move from one lattice site to another) in Sec. IH, Sec. IV 
presents an adaptation of Dirac's analysia23 to the case of 
a discrete lattice. The outcome is a gauge-invariant for- 
mulation of electromagnetism in empirical tight-binding 
theory. 

Although the theory derived in this way has many sim- 
ilarities with conventional tight-binding theory, there are 
significant differences as well. Not-all tight-binding mod- 
els can be made gauge invariant ;E2I this is possible only 
if the basis can be constructed from symmetrized co- 
ordinate eigenkets. In addition, local gauge symmetry 
imposes strong restrictions on the Hamiltonian matrix, 
which have the effect of sharply reducing the number of 
allowed Hamiltonian fitting parameters. Finally, unlike 
previous empirical tight-binding theories, the present ap- 
proach provides an explicit (discrete) wave function for 
the electron. 

The formalism derived here is applied to two semi- 
conductors with the diamond structure (Ge and Si) in 
Sec. For these systems, a basis of 15 orbitals per 
atom is shown to provide a satisfactory fit to the valence 
bands and the lowest conduction bands (up to about 
5 eV above the valence-band maximum). These results 
are comparable to those obtained from a 10-orbital basis 
proposed recently in ReL ^ The basis used here is 50% 
larger, but their modelEj cannot be made gauge invari- 
ant if intra-atomic coupling is included. Thus, it appears 
that some tradeoffs are necessary if gauge invariance is 
to be achieved. 
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II. COORDINATE MATRICES AND THE 
COORDINATE REPRESENTATION 



As mentioned above, the in tr a L-atnmip tj oi i pj ing used in 
existing tight-binding m©*ielJ£aOMBElGaO leads to a 
lack of gauge invariance.Ea This may be seen from a sim- 
ple sp^ model for a single atom. In this case, we know 
from atomic physics that there are coordinate matrix el- 
ements coupling the s and p orbitals: 

{s\x\p^) = {s\y\py) = {s\z\p^) EE c, (2.1) 

where c is real. In the basis {|s), \px), \Py), \Pz)}, the ma- 
trices representing x and y are therefore 
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(2.2) 



But this implies that x and y do not commute: 



xy — yx 



^0. 



(2.3) 



This means that the coordinate representation (consist- 
ing of simultaneous eigenkets of x, y, and z) does not 
exist. Even more important, it means that the theory 
cannot be gauge invariant. In a gauge transformation, 
the vector and scalar potentials transform as 



A^A + VA, 
and the state ket lib) transforms as 



IdA , , 



U = exp 



3A(x,t) 
he 



(2.5) 



where A is an arbitrary function of x = {x, y, z) and t. 
If a theory is gauge invariant, all physically measurable 
quantities must be independent of such transformations. 
But the expectation value (x) is a measurable quantity, 
and under a gauge transformation one has 



(U^xU), 



(2.6) 



where WxU 7^ x if A depends on y or z. Hence, no theory 
can be gauge invariant if a;, y, and z do not commute. 

Is there any way of achieving gauge invariance without 
setting c — 0? Perhaps the sp'^ basis is too small, and the 
situation might be improved by including more orbitals 
(d, /, . . . ). However, one soon finds that for any finite 
LCAO basis, the lack of gauge invariance persists. This 
follows directly from the Wigner-Eckart theorem — since 
x is a vector operator, it couples states with angular mo- 
mentum I to those with Z±l. Hence, any finite truncation 
of the basis results in non-commuting coordinates. 



Another possibility is to keep the same sp^ basis, but 
modify the coordinate matrix. The physical justification 
for doing so is the fact that the orbitals used in empiri- 
cal tight-binding theory are not-atomic orbitals; they are 
orthogonalized atomic orbitals.Ej Therefore, they do not 
have the full rotational symmetry of atomic orbitals; they 
have only the site symmetry of the crystal structure. For 
example, iiie atoms in a diamond crystal have site sym- 
metry TdSj Therefore, the orbital that was denoted \pz) 
above should really be written as jEfg), since it belongs 
to the Fis representation of TdS3 

However, the d orbital \dxy) also transforms as Irf 



Thus, in the Td group, the matrix element h = (Ffg \y\T 
is allowed, and the coordinate matrices (p^) become 
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0. 



(2.7) 



This yields 



xy — yx 






62_c2 











(2., 



which is equal to zero if 5 = ±c. 

Setting b — c would imply that the p and d orbitals 
have equal weight in the F15 states. This-is not as absurd 
as it sounds; Boguslawski and GorczycaEj have shown us- 
ing first-principles pseudopotential calculations that for 
the F15 states at the top of the valence band in GaAs, 
the probability of finding an electron in a cation d orbital 
is greater than that of finding it in a cation pxirbital. (In 
AlAs, the probability ratio is greater than 2.Efl) Thus, it 
is not unreasonable to assume that b and c have compa- 
rable magnitudes. 

If one sets b = c, then the coordinate operators have 
simultaneous eigenkets \x' ,y' , z'), which are given by 



\ccc) 
\ccc) 
\ccc) 
Iccc) 



- 5(|ri> 
= ^(iri) 

- ^(|ri> 
= 5(|ri> 



+ 



1^5) 



irfs)), 
|rf5», 
irfs)), 
irfs)). 



(2.9) 



Note that the coordinate eigenvalues are located at the 
corners of a tetrahedron. In fact, the linear combina- 
tions given in (|]^) are identical to the hybrid-bond 
orbitals used in analytical tight-binding theories,Lra al- 
though these are not ordinarily interpreted as exact co- 
ordinate eigenkets because the F15 states are assumed to 
be pure p orbitals. 

The procedure outlined above is rather clumsy; one 
simply modifies the coordinate matrices by trial and er- 
ror in an attempt to make them commute. One cannot 
predict in advance whether the attempt will succeed, and 



4 



in ge nera l it will not. However, the unitary transforma- 
tion (2.9) may be inverted to obtain 



iri) 



i(|ccc) 
i(|ccc) 
i(|ccc) 
i(|ccc) 



+ \ccc) 


+ ccc) 


+ ccc 


+ ccc) 


— |ccc) 


- |ccc 


— \ccc) 


+ ccc) 


- |ccc 


— Iccc) 


— Iccc) 


+ ccc 



which immediately suggests a more frui tful approach. 
The linear combinations given in Eq. ( 2.10| ) are just what 
one would obtain by starting with a single coordinate 
eigenket (say |ccc)) and using the symmetry operations 
of the t ptra hedral group to construct "symmetrized" 
orbitalstllcS that trarisform according to the irreducible 
representations of Tdx3 

Thus, in this alternative approach, the coordinate ba- 
sis is taken as fundamental, and the tight-binding ba- 
sis is merely a secondary alternative that is useful for 
reasons of symmetry. Since the existence of a coordi- 
nate representation is necessary for gauge invariance, no 
tight-binding basis can be made gauge invariant if it can- 
not be represented in terms of symmetrized coordinate 
eigenkets.EJ Hence, this symmetrization procedure pro- 
vides us with all possible gauge-invariant tight-binding 
models. 



The orbitals in Eq. (2.10) are useful as a starting point, 
but they cannot be interpreted as atomic orbitals be- 
cause they do not have inversion symmetry. To obtain 
more "atomic-like" orbitals, one can apply the symmetry 
operations of the cubic group Oh to the basis ket |ccc), 
which yields the orbitals 



iri) 



V8 



(I ccc) 



-f ccc) -|- ccc) 



I ccc)). 



— \fxyz) 

^^(|ccc) + I ccc) -I- I ccc) -I- I ccc) 
V 8 



ccc) — ccc) 



ccc) 



I ccc)). 



1 

^ 71 



(I ccc) — I ccc) — I ccc) -I- I ccc) 
— |ccc) + |ccc) + |ccc) — |ccc)). 



- \d 



xyl 



(2.11) 



V8 



(I ccc) 



ccc) — ccc) 



ccc) — ccc) 



I ccc)). 



Here two different labels are used: the representations 
of Oft ,£3 and the conventional atomic orbital notation. 
Other orbitals not given here may be obtained from cyclic 
permutation s of x, y, and z. Note that the orbital 
irfg) in Eq. ( 2.10 ) is the same as -^{\Pz) + \dxy)), whereas 
the Td orbital jEi) is just ^(|s) + \Uyz)). 



In the basis (2.11), the nonzero coordinate matrix ele- 
ments are 



A^IPx) = {Px\y\dxy) = {dxy\z\fxyz) 



(2.12) 



plus others given by cyclic permutations. The selection 
rules for x are thus the same as those in a spherically sym- 
metric atom — although, in any real atom, the matrix ele- 
ments ( 2.12| ) would not be numerica lly eq ual. This equal- 
ity occurs because the basis kets ( 2.11 ) are degenerate 



eigenkets of the radial coordinate r — \fx^ -1- 



To 



break the numerical equality, one would need to use basis 
functions with a linear combination of different radii. 

The above procedure may, of course, be applied to co- 
ordinate eigenkets other than |ccc). Below is a list of the 
representations obtained by applying the symmetry op- 
erations of Oh to several different "generator" eigenkets: 



1 000) 
1 100) 
1111) 
1110) 



Ti + ri5 - 
Ti + ri2- 



r25' 



-r2' 
r25' 



(2.13) 



Explicit basis functions for these representations are 
given in Appendix 

With these results, one can now construct a gauge- 
invariant tight-binding model simply by putting a set (or 
more than one set) of these "orbitals" on each atom in 
a crystal or molecule. In such a model, the coordinate 
operators x, y, and z commute by construction. However, 
one is no longer permitted to choose orbitals arbitrarily. 
The choices are limited to taking all of the orbitals in a 
set or taking none of them. As an example, one cannot 
discard the / orbitals in the basis generated by |110) 
without destroying the gauge invariance of the theory. 

This approach yields a tight-binding model with or- 
thogonal orbitals. Another approach is to define a grid 
of coordinates, some points of which are not uniquely as- 
sociated with individual atoms. One may still construct 
symmetrized orbitals in this case, but the orbitals are not 
orthogonal. This makes the tight-binding approach more 
difficult; however, one can simplify the theory by choos- 
ing a Bravais lattice for the coordinate grid, in which 
case the model may be viewed as a discrete pseudopoten- 
tial model. Applications of both the pseudopotential and 
tight-binding approaches are considered below in Sec. 



III. TOPOLOGY OF THE LATTICE 

As we have seen, the most general gauge-invariant 
tight-binding basis consists of a set of discrete coordinate 
eigenkets, which will be referred to as a lattice. Such a 
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lattice is generally not periodic. In order to apply the 
principle of local gauge symmetry to such a system, one 
must be able to calculate the change in phase thjat, occurs 
along any specified path in coordinate space.E3'E3 Thus, 
the first step is to define what is meant by a "path" in a 
discrete coordinate system. 

In general, a path is just an ordered sequence of points. 
In a continuous coordinate system, neighboring points in 
the sequence must be separated by an infinitesimal dis- 
tance. This defines the topology of the system, in which 
points are linked together only if they are adjacent in co- 
ordinate space. It is desirable to define the topology of 
the discrete lattice in a similar way. 

One way to do this is to construct-a Voronoi polyhe- 
dron around each site in the lattice.cZl A Voronoi polit= 
hedron is just the region in space closest to that point;E3 
if the lattice is a Bravais lattice, the polyhedron is the 
same as a Wigner-Seitz cell. In mathematical terms, the 
Voronoi polyhedron fli for site is the set of points x 
such that |x — Xi| < |x — Xj| for all j ^ i. The topology is 
then defined by the following rule: If Hi and share a 
surface with area Sij > 0, the sites x^ and x^ are linked 
together; otherwise, they are not. 

In some cases, it may happen that two Voronoi poly- 
hedra share only a point or a line, in which case Sij = 0. 
The linking algorithm presented in Ref. ^ does not con- 
sider this possibility (because Ref. ^ deals only with ran- 
dom lattices, for which the probability of such an event 
is zero). For certain applications, it is useful to include 
links between such sites,c3 but we shall see below that 
these links should be excluded in the present situation. 
Thus, only adjacent sites whose Voronoi polyhedra share 
a surface with Sij > are linked. 

A path in the discrete lattice is then just an ordered 
sequence of linked points, and a closed path is one whose 
first and last points are the same. By definition, every 
edge of a Voronoi polyhedron is equidistant from three or 
more lattice sites, all of which lie in a plane perpendicular 
to the given edge. These sites are closer to this edge than 
any other sites. The links between these sites form a 
closed path, and the area bordered by the links is called 
a "plaquette." There is a one-to-one relationship between 
the plaquettes and the edges of the Voronoi polyhedra. 

The plaquettes partition all of coordinate space into 
nonoverlapping volumes called (Delaunay) cells. Each 
cell is uniquely associated with one corner of a Voronoi 
polyhedron. The partitioning of space into cells is re- 
ferred to as a Delaunay tessellationHa 

A general algorithm for calculating the geometry of 
Voronoi polyhedra, links, plaquettes, and cells is pre- 
sented in Appendix ^ The expressions derived there 
will be of use in what follows. 



IV. LOCAL GAUGE SYMMETRY ON AN 
ARBITRARY DISCRETE LATTICE 

Christ, Friedberg, and LeeS'H'Ei have developed a 



theory of local gauge symmetry on a random lattice. 
This section presents a modified version of their theory, 
with special emphasis on the implications of the prin- 
ciple of local gauge symmetry for tight-bindipg theory. 
The presentation follows Dirac's approachJ2j in which 
the existence of electromagnetic fields is "derived" as a 
straightforward consequence of a degree of freedom (non- 
integrable phases) possessed by any quantum-mechanical 
system that can be represented in a coordinate basis. 

A. Electromagnetism is a nonintegrable phase 

In a discrete coordinate basis, any ket vector may be 
expressed as 

|V') = ^c.|x,), (4.1) 

i 

where |xi) is a coordinate eigenket, which is assumed to 
be normalized such that 

(x,|x,)=<5,,. (4.2) 

Dirac's starting pointS is the fact that physical predic- 
tions in quantum mechanics are ultimately expressed in 
terms of probabilities of the form Kf/jj^/;')!^, where the 
probability amplitude {ipl'ip') is given by 

(V'l^')-Ec*^^ (4-3) 

i 

The probability is obviously well defined even when the 
overall phase of {ip) has no definite value. This degree of 
freedom is referred to as global gauge symmetry. 

The existence of global gauge symmetry raises the 
question of whether it is necessary for the local proba- 
bility amplitude q to have a definite phase. In other 
words, suppose we write 

Q = 5,e'ft, (4.4) 

where the phase of bi is well defined (to within an integer 
multiple of 2tt), but Pi is a nonintegrable function — that 
is, the change in Pi around a closed path can take on 
any value. In this case one can see that | ■;/)') p is well 
defined only if the change in Pi around any closed path is 
the same for all states \^) and \ip') (to within an integer 
multiple of 2tt, which is absorbed into the definition of 
bi). But anything that is the same for all states can be 
viewed as a physically real part of the dynamical system. 
Since the present system consists only of a single point 
particle, these nonintegrable phases must represent a field 
of force acting on the particle. 

The principle of local gauge symmelcy is therefore de- 
fined by the following two postulates:Lj (i) The physical 
predictions of the theory must be unambiguous, (ii) The 
phase of Ci at any point in space and time need not be 
well defined; only the change in phase between linked 
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points must be definite. As shown above, these postu- 
lates entail that the change in Pi around any closed path 
must be the same for all states. According to postulate 
(ii) , this change is fixed for any path by the change in (3i 
between two linked points in space: 



= A - /3, 



(i linked to j), 



and in time: 



(4.5) 



(4.6) 



Since (3i is nonintegrable, Kij and are independent vari- 
ables. These two quantities are the fundamental dynam- 
ical variables that arise from the principle of local gauge 
symmetry. It will now be shown that Hij and can be 
interpreted as potentials for the electromagnetic field. 

One possible closed path involves a space displacement 
dij = Xi — Xj and an infinitesimal time displacement di, 
followed by dji and —dt. The change in phase around 
this path can be used to define (tentatively) an electric 
field variable 



Eji - - 



Aj — A^ 



dji 



(4.7a) 



where c?,- 



If the index £ is used to label the links, 



one may write this in the simpler form 



h A\g — ki 
e di 



(4.7b) 



Here Ee is interpreted as the component of the electric 
field in the direction = dji] the components perpen- 



dicular to dg are not defined. Equation (4.7) takes a 
familiar form if expressed in terms of the potentials 



since then 



— —^n 

e 



Ef 



he Kl 

Ai = 

e dp 



c dt 



(4.8) 



(4.9) 



Here the notation dAg/dt indicates that di is to be held 
constant during the differentiation. 

Another type of closed path is an elementary plaquette 
q constructed from links t in coordinate space (see Sec. 
Ill and Appendix The change in phase around the 
perimeter of the plaquette may be used to define the 
magnetic field 



(4.10) 



where Sq is the area of the plaquette [see Eq. (B24)]. 
Summing Eq. (4.10) over the (closed) surfaccref a cell c 
leads immediately to the "no monopoles" lawE3 



(4.11) 



Likewise, summing E^di around the perimeter of a pla- 
quette gives Faraday's law: 



Edt 



1 d 

cdt 



(BqSq). 



(4.12) 



The other two Maxwell equations can be obtained from 
the Lagrangian L = Lf + L^, where Lf is the electromag- 
netic field Lagrangian 



Lf = ^y3E'i^i~^y 



(4.13) 



Here ilf, — ^Sed^ 



is the volume of link £, where Si = Sij 
is the area of the surface shared by the Voronoi polyhedra 
for sites i and j (see Appendix |b|). Likewise, V,q — \Sqdq 
is the volume of plaquette q, where dq is the length of 
the Voronoi polyhedron edge corresponding to q. Equa- 
tion ( 4.13| ) is, just a discrete version of the standard field 
LagrangiarH JiE^ — B^)d'^x; the only apparent dif- 
ference is an extra factor of 3. This factor cancels the 
factor of i in the definition oiilp 



and r2„, thus leading to 



t he co rrect Maxwell equations below. It appears in Eq. 
(4.13) because the standard Lagrangian is expressed in 
terms of E'^ = E^ + Ey + E"^, whereas Ef includes only 
the component of E in the direction of d^. 

The electronic term in the Lagrangian, which includes 
the field-particle coupling, is 

Le = ^*^^'^ (4.14a) 

= T E(^*^^ - ^*^*) - E (4-14b) 



Here Hij = 



(x, |i7|x.,) is the Hamiltonian in the absence 



of electromagnetic fields, while 



Hij — Elij i 



^^^^^^ 



(4.15) 



The first expression ( 4.14a ) for Le has exactly the same 
form as the Lagrangian in the case of no electromagnetic 
fields. This expresses the fundamental physical content 
of the principle of local gauge symmetry — that the in- 
fiuence of the field upon the particle can be expressed 
entirely in terms of the nonintegrable phase of the prob- 
ability amplitude Ci = bie^^\ In the second expression 



( 4.14b| ) for Le, all of the nonintegrable phases are col- 



lected together in the eflECctive Hamiltonian Hij. This is 
the usual approach, in which the probability amplitude 
bi has a well-defined phase, and the field appears only in 
the Hamiltonian. 



The Hamiltonian (4.15) appearing in the Lagrangian 
( 4.14b| ) depends upon the phase difference /3i — f3j . This 
phase difference is not well defined unless the sites i and 
j are linked. But according to postulate (i) above, all 
physical predictions of the theory must be unambiguous. 
Hence, the principle of local gauge symmetry demands 
that 



Hi, 







{i not linked to j) 



(4.16) 
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Local gauge symmetry therefore imposes constraints not 
found in conventional tight-binding models. 

Note that the Lagrangian L is gauge invariant by con- 
struction. In other words, both Lf and Lg are invariant 
under the gauge (phase) transformation 



h 
A, 



+ Xi 



(4.17) 



^ Xi Xj 1 



where Xi is an arbitrary integrable function. 

Given the above Lagrangian, the Euler-Lagrange equa- 
tion for Xi or (f>i is just Gauss's law 



(4.18) 



where qi — eb*bi is the charge on site i. The correspond- 
ing equation for Kg is the Ampere-Maxwell equation 



Bqdg — z~n(^^^^) ~ — -^^ 



d„eSe 



(4.19) 



where If 



'ji = {2e/ h)l ia{b* H jibj) is the current from 
site i to site j. Summing ( [1.1S| ) over all links tha|_contain 
a given site i yields the charge conservation lawEj 



= 0. 



(4.20) 



Thus, we see that Xi and Kg can be given a consistent in- 
terpretation as discrete electromagnetic potentials, since 
the above equations are in full agreement with macro- 
scopic (i.e., long- wavelength) electromagnetism. 

In some applications of Voronoi polyhedra, it is useful 
to link sites i and j whose polyhedra share only a line 
or point, hence St = Sij = 0.E3 For such links, the link 
volume Vlf^ = ^Sidi is zero, so the electric field Ei does 
not contribute to the field Lagrangian ( 4.13|), G auss's 
law ( 1.1^ ), or the Ampere-Max well eq uation ( 4.19| ). The 
magnetic- field contribution to ( 4.19 ) likewise vanishes, 
because Si = Q. The current through such a link must 
therefore be zero, which can only be true in general if the 
Hamiltonian matrix element Hij vanishes. Links with 
Sij — are consequently devoid of any physical signifi- 
cance, and there is no loss of generality if one excludes 
them at the outset by linking only sites with Sij > 0. 

Returning to the Lagrangian L, the Euler-Lagrange 
equation for b* is just the Schrodinger equation 



ihbi = ^ Hijbj 



(4.21) 



Note the strong simi lari ty between this result and the 
Peierls substitution (|l.4|). The main difference is that 



(4.22) gives the Hamiltonian in the coordinate represen- 
tation, not the tight-binding representation. 

If \Kij\ ^ 1 (i.e., if the field is weak or the lattice 



spacing is small), the Hamiltonian (4.22) reduces to 
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Hij ~ Hij — — • pij 



^A,,+eq^A,. (4.23) 



Here a vector potential has been defined as Aij — Aij dij , 
while the momentum operator is given by 



Py = ^dyJ/y = — (xi - 



in ih 

which is the s ame as the kinematic mo mentu m defined 
above in Eq. (p|). The Hamiltonian ( 4.22 ) therefore 
clearly gives the correct first-order A • p coupling. We 
shall see below that the dimensionless quantity 

1 Tfl 

\j = -^d,, • p„ = -j^d^H^j (4.25) 

can be viewed as a geometric weight factor that gives the 
correct coupling also. 



B. Geometric definition of momentum 
and kinetic energy 

Up to this point, little has been said about the struc- 
ture of the Hamiltonian Hij. Within the bounds of the 



restriction (4.16), Hij may be treated as an arbitrary fit- 
ting parameter. However, in some circumstances it may 
be desirable to reduce the number of fitting parameters 
by using a theoretical formula for Hij that would repro- 
duce the Schrodinger equation in the limit of zero lattice 
spacing. 

Let us start by considering the momentum operator, 
which will be defined in this section as the canonical mo- 
mentum p = —ihS7. A discrete expression for the gradi- 
ent operator may,be obtained from the integral definition 
of the gradients 



V/(x) = lim 



(4.26) 



where dS is a surface element pointing outward from fJ. 
Now the limiting volume in a discrete lattice is the vol- 
ume r2i of the Voronoi polyhedron for site x^. On the 
surface Sij shared by il, and ilj, the value of /(x) may 
be taken to be ^[/(xi) -f /(xj)]. Hence, the discrete gra- 
dient may be defined as 



Since bi is an ordinary probabihty amplitude with a V/(xi) = — ^ -[/(x^) + /(xj)]Sji, (4.27) 

- - - - ~ - \li ' 



well-defined phase, Hij must be the Hamiltonian in the 
prese nce of electromagn etic fi elds. With the restriction 
( 4.16 ), one can express ( 4.15 ) as 



where Sj^ — Sjidji. Now since 



Hi, 



H,j exp{-iKij) + hXAj 

Hij e'xplieAijdij / he) + eipiSij. 



(4.22) 



dS = 0, 



(4.28) 



8 



the term involving /(x^) drops out, leaving only 

V/(x,)--J-^/(x,)S,,. (4.29) 



An alternative derivation of this result is given in Eq. (17) 
of Ref. H. 

The caijxinical momentum operator p may therefore be 
defined aM 



{^^^\P\'P) 



(4.30a) 
(4.30b) 



— 1/2 

Here the basis kets |Xj) = fl.^ \xi) are chosen to satisfy 
"(5-function" normalization 

(X.|X,) 



(4.31) 



in contrast to the usual kets Ixi), which are normalized 



to unity [see Eq. (4.2)]. The normalization (4.31) is used 
here because it agrees (in the limit Oj — » 0) with the 
i5-function normalization of cont inuous coordinate, eigen- 
kets, upon which the definition ( 4.30a ) is basedH 
Substituting \ip) ~ \X.j) in Eq. ( 4.30b ) then gives 



(x.IpIx,) 



ihSij 

2 ^ J 



which clearly satisfies 

5](X.|p|X,)l], =0. 



(4.32) 



(4.33) 



Replacing |Xj) = il^ ' |xi) in ( ^4.32[ ) then yields the de- 
sired result 



(X,;|p|Xj 



ihSi 



2 i j 



(4.34) 



Note that this matrix is Hermitian, because 8,^ = — 



If the kets jx^) are used in Eq. (4.30a) above, a non- 
Hermitian canonical momentum is obtained. 

There is some question as to whether this definition of 
p should be referred to as "canonical," because it does 
not satisfy the canonical commutation relations. In a 
continuous coordinate basis, the canonical momentum 
satisfies 



(x'|[a;",/]|x") =zM„05(x'-x"), 



(4.35) 



where a and /3 are Cartesian components of the given 
vectors. The corresponding equation in the discrete basis 
is 



(4.36) 



which obviously does not agree. Note, however, that 



^(x,|K,/]|x,)f7,f7, = yE'^S-^i 



(4.37) 



where fl is the total volume, and the second equality 
is proved in Eq. (11) of Ref. ^ This agrees with the 
relation 



(x'|[a;",/]|x")dVdV' = ihSa^pQ (4.38) 



in the continuous basis. Hence, Eq. (4.37) is as close as 
one can come to a carianical commutation relation in a 
general discrete basis .EjS 

A similar definition may be used for the kinetic energy 
operator — ?i^V^/2m. The integral definition of the 
divergence)^ 



V • F(x) = hm 



gives the Laplacian 

VV(x) = lim 
the discrete form of which is 



(4.39) 



(4.40) 



(4.41) 



An alternative derivation of this result is given in E q. (1 2) 
of Ref. E^. The procedure used above in Eqs. (4.3C)- 



(4.34) then yields the kinetic energy operator 



(xj|T|xj) 



Si 



2m dij 



Sik 

k 



which satisfies [cf. Eq. ( |4.33| )] 

^(X,|T|X,)r!, 



(4.42) 



(4.43) 



Note that for i ^ j, (xi|T|xj) decreases continuously to 



zero when S*,. 



0. This ensures that the Hamiltonian 



is a continuous function of the lattice coordinates, even 
as new links are formed and old ones are broken. 

Such continuity is also desirable when the Hamilto- 
nian is determined empirically—especially for applications 
(such as molecular dynamicsBo) in which the atomic po- 
sitions vary with time. This can be achieved by defining 
the nonlocal elements of the empirical Hamiltonian as 



(x,:|H|x,) = 



Sj' 



2m d-ij^^irij 



{i^j), (4.44) 
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where the fitting parameter fij is a continuous, nonsin- 
gular function of the lattice coordinates. 

Note that the operators p and T do not satisfy T = 
p^/2m, because p^, unhke p and T, couples sites that 
are not linked. However, p and T are related by 



m 

Pij = —dijTij 



(4.45a) 



C. Spin 

The theory presented thus far has been for a particle 
with spin zero. Particles with spin i may be described 
using a discrete version of the Dirac Hamiltonian for a 
free particle: 



H ~ CO. ■ p + (3nic , 



(4.51) 



in 



(4.45b) 



Thus, any Hamiltonian of the form H = T + V{x.), where 
y is a local potential, satisfies Eq. (1.1). Hence, for such 
a Hamiltonian, the canonical momentum p = —ihV used 
in this section agrees with the kinematic momentum de- 
fined earlier. 

Now let us examine the dimensionless factor A,-, de- 



fined above in Eq. ( [4.25| ). U H = T +V, this becomes 

3 j n 



i-^ij — 



(4.46) 



where ftij = \Sijdij is the volume of the link between 
sites i and j (see Appendix The factor Ay appears 
in the term in the Hamiltonian ( |4.23|) , which will be 
referred to as H2. In a continuous coordinate basis, H2 
is given by 



2mc^ 

which means that it satisfies 
[■x!\H2W)A\'Ay = 



A2(x')(5(x' -x"), 



(4.47) 



A2(x')dV. (4.48) 



The corresponding equations for the discrete basis are 

J2 



(XilFzlX,-) = ^ — ^ 



(4.49) 



and 



2mc^ 



(4.50) 



The second equality in ( 4.5Cl| ) was obtained by noting 
that a sum over i and j covers each link i twice. The 



only apparent difference between Eqs. (4.48) and (4.5C) 
is a factor of 3. This appe ars fo r the same reaso n tha t 
it does in the Lagrangian (4.13) — i.e., i n Eg. ( 4.4§| ) 
refers to + Ay + A^, whereas Aj in Eq. ( 4.50 ) refers 
only to the compo nent o f A in the direction of d^. 

Therefore, Eqs. ( 4.48 ) and ( 4.50 ) are the same in the 
limit of zero lattice spacing, and the factor Aij is simply 
a geometric weight factor th at pr ovides the correct A^ 
coupling in the Hamiltonian (4.23). 



where ct and /3 are Dirac's 4x4 matrices. The momentum 
operator p can either be calculated from geometry or 
fitted to experiment. In the presence of electromagnetic 
fields, the Hamiltonian becomes 



H = CO. ■ TV + Pmc^ 



where [cf. Eq. ( [4.22| )] 

T^ij = Pij exp(-iKij). 



A nonrelativistic Hamiltonian may be obtgiuae 
plying a Foldy-Wouthuysen transformationE3E 
( [4.52 ), which yields 



(4.52) 



(4.53) 

by ap- 
to Eq. 



1 



(cr-7r)* + e0 (4.54) 
cr • TT, {[cr ■ TT, ecj)] + ihcr ■ tt)]. 



where cr is the Pauli spin matrix, and all terms of or- 
der (v/c)* have been included. This Hamiltonian cou- 
ples sites that are not linked, but there is no ambiguity 
because the Dirac equation is taken as fundamental. 

If we assume for simplicity that the lattice coordinates 
do not depend on time, then 



where 



V^J = e( 



incr ■ TTiii = cr ■ tt. 



ij Vij , 



(4.55) 



(4.56) 



is the difference in potential energy of sit es i and j due 
to the electric field. The last term in Eq. ( 4.54 ) therefore 
consists of the Darwin term 



1 



■ {Vk^ + Vkj)n,k ■ TTfej (4.57) 



plus the spin-orbit coupling 



F{^° — —— 



"^iVkt + Vkj)cr ■ {TVik X TTfej), (4.58) 



where the identity 

(cr • a)((T • b) = a • b + icr • (a X b) 



(4.59) 



has been used. Now the main contribution to spin-orbit 
coupling comes from the atomic cores, where the poten- 
tial energy and wave function vary rapidly. However, in 
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any basis of reasonable size, the lattice imposes a wave- 
length cutoff that eliminates such rapid variations. The 
potential 0; must therefore be viewed as a pseudopoten- 
tial, not a true atomic potential. Hence, for p ractical 
purposes, in the spin-orbit Hamiltonian (4.58) should 



be treated as a fitting parameter that is ind ependent of 
the value used for the first ecf) term in (4.54). 

The first two terms in the Hamiltonian (4.54) are 



kinetic-energy terms, which may be rewritten using 

(cr • 7r)jj- = ^[TTifc • TTfcj + i(T ■ (TTjfc X TTfcj)], (4.60) 



in which the second term describes the intrinsic magnetic 
dipole moment of the particle. For a general lattice, this 
term is not zero even when there is no electromagnetic 
field, because different components of the momentum op- 
erator do not commute (i.e., pxp ^ 0). This follows from 
the fact that there is generally no more than one site k 
linked to both i and j, and for that i and j, Tpik x Pfcj is 
generally not zero. 

However, if the lattice is a Bravais lattice, then pxp 
is always zero. This follows from the fact that every site 
in a Bravais lattice is identical, so for a given nonzero 
Pife X Pfej, there is always another site I with p,; = p^j 
and pij = pifc, hence Pifc x pfc^- -I- p^; x pij = 0. Clearly 
one also has dik x d^j + du x dy = 0, so the sites i, k, j, 
and I lie in a single plane. If i is not linked to j and k is 
not linked to Z, then i, fc, j, and I form a single plaquette 
q, which has the shape of a rectangle. Otherwise, they 
form two triangular plaquettes. 



If the momentum operator is given by Eq. ( 4.34 ), then 
t he in trinsic magnetic dipole term in the Hamiltonian 
( 4.54 ) for such a Bravais lattice is 



rrinag 



eh j SikSkjSq 



8mc \ dikdkj^i 



(4.61) 



where the weak- field approximation \Kij\ <C 1 has been 
used, and the direction of is that of S^. If the sites i, 
k, j, and I form a single rectangular plaquette, then Sq is 
the area of that plaquette; otherwise, it is the combined 
area of the two triangles (in which case is the average 
magnetic field of the two plaquettes). 

As an example, consider a simple cubic lattice with 



lattice constant 



for which Sik — Shi = S„ — 



dik = dkj = a, and Q,,, = 



In this case, the factor in 



parentheses in Eq. ( |4.6l| ) is unity, and H^^^ couples sites 



on opposite corners of each plaquette (with dij = ^/2a). 
By comparison, the dipole term in the continuum Hamil- 
tonian is given by 

-^[<T . (p - eA/c)r = -^(p - eA/c)2 - ^cr ■ B. 

(4.62) 

The numerical factor in front of thi s dip ole coupling is 
four times larger than that in Eq. (4.61). This occurs 
because (4.61) couples each site i to four other sites j. 



V. APPLICATION TO TETRAHEDRAL 
SEMICONDUCTORS 

This section considers several different methods of im- 

Spin is ne- 



plementing the theory developed in Sec. [V 



glected in all of the applications that follow. 



A. Discrete pseudopotential method 

The simplest geometry occurs when the lattice sites 
Xi are chosen to lie on a Bravais lattice. One possible 
appro ach in this case is to use the geometric expression 
(4.42) for the kinetic energy T, and assume that the po- 
tential energy V is local. This approach will be referred 
to as the discrete pseudopotential method. 

If X lies on a Bravais lattice (the subscript i is omitted 
here) , one may define a reciprocal lattice as the set of all 
vectors g such that g • x = 27r x integer. The volume 
of a primitive cell in the direct lattice is denoted ujq, 
while that of a primitive cell in the reciprocal lattice is 
= (27r)3/wo. 

In a crystal, the Hamiltonian will be periodic with re- 
spect to some larger Bravais lattice whose sites are de- 
noted by R, where R G {x}. One may then define a 
corresponding reciprocal lattice as the set of all vectors 
G such that G • R = 27r x integer. The volume of a 
primitive cell in R space is fio = nLDQ, where n is an in- 
teger, and the volume of a primitive cell in G space is 
= (27r)3/f^o. 

Periodic boundary conditions may then be imple- 
mented over an even larger Bravais lattice whose sites 
are denoted by L, where L S {R}. The corresponding 
reciprocal lattice vectors are denoted k, where k • L = 
27r X integer. The volume over which periodic boundary 
conditions are applied is SI = -/VSIqj where N is an in- 
teger, and the volume of a primitive cell in k space is 
SI* — {2tiY /VL. Note that according to the above defini- 
tions, G e {k} and g € {G}. 

Coordinate eigenkets in this system are denoted |x), 
where |x -I- L) = |x) due to the periodic boundary con- 
ditions. The orthogonality and closure relations in this 
basis are therefore 



L 



(5.1) 
(5.2) 



In a system with periodic boundary conditions, the co- 
ordinate operator is not well defined; only periodic func- 
tions of the coordinate are permitted. Therefore, the 
definition of the kinematic momentum must be slightly 
modified. Instead of (1.1), one has 



(xlplx') = — (x - x')(x|i/|x') if x e and x' G 

ih 



(5.3) 
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with (x|p|x') = (x + L|p|x') otherwise. 

Another useful representation is the crystal momentum 
representation 



e^-'^lx) 



(5.4) 



where Af = nN = ^I/luq. The corresponding orthogonal- 
ity and closure relations are 



(k|k'>=^5k-k',g, 
g 

E Ik>(k| = l. 



(5.5) 
(5.6) 



In the crystal momentum representation, a periodic 
Hamiltonian (x|if|x') = (x + R|i/|x' + R) couples only 
those states that differ by a reciprocal lattice vector G: 



(k'|JT|k) = E -^k-.k+cik + G|iJ|k), 



(5.7) 



G 



where 



(k + G|i7|k) = E e-*(''+°)-"'(x'|i7|x)e''"". 

(5.8) 



The kinematic momentum (5.3) also satisfies Eq. ( p.7[ ). 
Its matrix elements that are related to those of H by 



(k + G|p|k) = -Vk(k + G|i7|k) 



(5.9) 



On a Bravais latti ce, th e kinetic energy ( 4.42 ) and 
canonical momentum (4.34) are translationally invariant: 



(x|r|x')=r(x-x'). 



(5.10) 



where r(x + L) = T(x). The matrix elements of T are 
therefore given by 

(k|T|k')-r(k)E'5k-k',g, (5.11) 



where 

T(k) = J2 7'(x)e"*-^. (5.12) 

For tetrahedral semiconductors with the diamond or zinc- 
blende structure, it is convenient to use a cubic lattice for 
the grid {x}. Expressions for the link distances dij and 
surface areas Sij are given in Appendix |C| for the sim- 
ple cubic, body-centered cubic, and face-centered cubic 
lattices. The resulting kinetic-energy operators given by 



Eqs. ( p^ and ( |5l^ ) are 
2ft2 



rsc(k) 

7bcc(k) 
rfcc(k) 



ma 



2 |^sin^(^fc2;a) + sin^(^fcj,a) + sin^(^fc2a)]. 



2 2 {^[-'^ ~ cos(ifc^a) cos(iA;j^a) cos{^kza)] 
+ [sin^(ifca;a) + sin^i^kya) + sin^^fc^a)]}, 

^[3 — cos(iA;„a) cos( ikza) 

— cos{^kza) cos{^kxa) 

— cos(^A:j.a) cos(^fcj^a)], (5.13) 



all of which reduce to r(k) ~ h?k^ /2m when fca ^ 1. 
Here a is the lattice constant of the grid {x}, which is 
some integer fraction of the lattice constant oq of the 
crystal lattice {R}. 

A canonical momentum operator p(k) corresponding 
to Eq. (4.34) may be defined in a similar manner. This 
operator is given by 



p(k) - -VkT(k), 



(5.14) 



which follows from Eq. (4.45). Note that this resu lt is 
just a special case of the kinematic momentum (5.9). 

The matrix elements of a local pe riodic potential 
y(x) = F(x + R) are given by Eq. (p/7|), where 



(k + G|^|k) = i E 



iG-. 



(5.15) 



is independent of k. It will be assumed here that is a 
superposition of local atomic pscudopotentials: 

where w^(x) = t;^(x-|-L) is the pseudopotential for atom 
/i, whose position in the unit cell Q.q is given by r^. In 
this case 



(k + G|l/|k) 



rE-M(G). 



Na 



(5.17) 



where v^{G) is the atomic form factor 



Na 



E ^/^W"^ 



iG-: 



(5.18) 



xGH 



and Na is the number of atoms in the unit cell SIq- 

The main practical difficulty in implementing the dis- 
crete pseudopotential method is that T(k) is not a good 
approximation to the continuum kinetic energy 



Tco„t(k) 



21,2 



2m 



(5.19) 
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FIG. 1: Energy band structure of GaAs calculated by the 
discrete local pseudopotential method using an fee gri d wit h 
a = |ao. Solid lines: fee kinetic energy from Eg. (3.13). 
Dotted lines: continuum kinetic energy from Eq. 



a constant factor 7r2(2 + \/2)/32 = 1.053 so that rfcc(k) 
matches Tcont(k) at G (111) and G = (200). The 
pseudopotential form factors ( |5.18| ) for this calculation 
were taken from Ref. ^ No attempt was made to fit the 
energy bands by modifying the form factors; the purpose 
of this figure is merely to demonstrate the close similarity 
between the discrete fee band structure and the contin- 
uum band structure. Slight adjustments in the model 
parameters would likely give an even better agreement. 

The main problem with this result is that a = ^ao cor- 
responds to a basis size of 512 grid points per primitive 
unit cell fJo- This is unattractive in comparison to the 
basis dimensions of approximately 100 plane waves that 
are typically used in empirical pseudopotential calcula- 
tions for tetrahedral semiconductors. However, changing 
the fee grid to a = |;ao (i.e., 64 grid points per primitive 
cell) makes it impossible to achieve a satisfactory fit to 
the band structure using local pseudopotentials. A good 
fit is possible only if the nonlocal Hamiltonian matrix 
elements are treated as fitting parameters. But in that 
case, one can reduce the basis dimensions even further 
by using the tight-binding approach. 



unless fca ^ 1. This means that the lattice constant a 
of the grid {x} must be significantly smaller than the 
lattice constant oq of the crystal lattice {R}. To obtain 
one grid point at each atom in the diamond structure, 
a must satisfy a = ao/21 (for a bcc grid) or a = a^/Al 
(for sc and fee), where I is a positive integer. Numerical 
accuracy generally requires Z > 1, as shown below. 

The shape of the Brillouin zone is also important. If 
the shape of the Wigner-Seitz cell for Wq is not con- 
gruent with the shape of the Wigner-Seitz cell for JIq, 
then T'(k) will deviate from Tcont(k) more rapidly in 
some directions than others. This can lead to significant 
qualitative errors in the kinetic energy. For example, 
in diamond, the lowest continuum eigenvalues at T are 
rco„t(G) = ^2^2/2™, where G = (000), (111), (200), 
and (220) (in units of 27r/ao). But for a sc grid with 
a = ^oo, the value of rsc((200)) = jmc? is actually 
lower than that of Tsc((lll)) = ZT? jma}. The correct 
ratio Tcont((200)) = |Tcont((lll)) is only approached in 
the limit of very small a/ao, making the sc grid a poor 
choice for diamond. 

The natural choice for diamond is the fee grid, since 
its Brillouin zone has the same shape as that of diamond. 
Indeed, one has Tfcc((200)) = |Tfcc((lll)) even for the 
maximum grid size a = \aQ. The only problem here is 
that Tfcc((220)) = |rfcc((200)), which is not sufficiently 
close to the correct ratio Tcont((220)) = 2rcont((200)) to 
make a = \aQ a satisfactory choice for numerical calcu- 
lations. The next possibility is a = ^oq, which yields 
Tfcc((220)) = (I + ^)Tfc,((200)) ~ 1.85Tfce((200)). 

The energy band structure for GaAs calculated using 
an fee grid with a = |ao is given in Fig. ^ The fee ki- 
netic energy obtained from Eq. (5.13) was multiplied by 



B. Tight-binding method 

In the tight-binding approach, the grid points are no 
longer restricted to lie on a Bravais lattice, and all of the 
Hamiltonian matrix elements are treated as fitting pa- 
rameters. In addition, it is assumed here that the model 
is constructed using the symmetrization procedure de- 
scribed in Sec. ||, so that a distinct set of orthogonal or- 
bitals is associated with each atom. The objective then 
is to find the smallest coordinate basis that provides a 
physically reasonable model of a given system. 

The basis kets in the tight-binding approach will be 
written as \a, R-t-r^^), where R is a lattice vector for the 
Bravais lattice over which the Hamiltonian is periodic, 
and is the position of atom ^ within the primitive 
unit cell Op. The se quantities are defined exactly as they 
were in Sec. V A the vectors L, G, and k are also defined 



in the same way. The label a may refer to a coordinate 
Xq within the atom, in which case 



\a,R + Tn 



- R + T„ 



(5.20) 



is just a coordinate eigenket. However, a may also be 
used as a symmetry label for an atomic orbi tal th at is a 
symmetrized linear combination of the kets ( |3.20 ): 



,R 



(5.21) 



In either case, the basis is orthogonal: 
(a,R + T,Ja',R' + v) = 5aa'<5^M' '^R-R-'^l. (5-22) 
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and complete: 

^ |a,R + r^)(a,R + r^| = 1. (5.23) 

a 11 Hen 

In pjeriodic systems, it is convenient to define the Bloch 
sumsEJ 

|a,/x,k) = ^ E e^'^-(^+"-V,R + r^), (5-24) 
Hen 

which are also orthogonal and complete: 

(a, ^, k|a', /i', k') = 5aa'5titi' ^ <5k-k',G e"*'^"^^' , 



G 



EE E |a,^,k)(a,//,k| = 1. 

a 11 keOp 



(5.25) 



(5.26) 



If the Hamiltonian H is invariant with respect to lattice 
translations R, then its matrix elements in the Bloch 
basis are 



(a, /i, k|_ff|a', //', k'^ 



(a,/i,k|ff|a',^',k) (5.27) 



where 



X e 



(5.28) 



The kinematic momentum operator (4.24) is related to 
this Hamiltonian by 

TTl 

(a,/Lt,k|p|a',/i',k) = — (x^ - Xc,/ + iVk) 
in 

x(a,/i,k|if|a',Ai',k),(5.29) 

where the original bas is was assumed to be given by 
(5.20). Note that ( ^.29| ) has the form of an intra-atomic 
matrix element (proportional to Xq — Xq,/ ) plus an inter- 
atomic matrix element (pr opo rtional to iVk). The zero- 
parameter model of Eq. (1.3) is obtained in the limit 
Xq 0. 

Let us now consider specific examples of H for tetrahe- 
dral semiconductors. The simplest tight-binding model 
for the diamond or zinc-blende structure consists of a 
single s orbital per atom, which is obtained by putting 
one c oordinate eigenket at each atomic position [see Eq. 
( 2.13 )]. In this model, as shown in Appendix |y, each 
atom is linked to 4 nearest neighbors and 12 second- 
nearest neighbors. More distant linkages do not exist 
because the Voronoi polyhedra for these atoms do not 
touch one another. 

Such a simple model is, of course, unable to describe 
even the qualitative features of tetrahedral semiconduc- 
tors. The simplest conventional tight-binding model that 



TABLE I; Number of free parameters for different tight- 
binding models in the diamond structure. The upper half 
of the table refers to conventional tight-binding models, while 
the lower half refers to models obtained from symmetrized 
coordinate eigenkets. 



Model 






Parameters 




Basis 


Size 


On-site 


1 MM 




3NN 




4 


2 


4 


7 


7 


sp^s* 


5 


4 


7 


11 


11 


sp^d^ 


6 


3 


7 


13 


13 




10 


7 


17 


33 


33 


|iii> (r^) 


4 


2 


1 


1 





1100) 


6 


2 


1 


1 





illl) (Oh) 


8 


3 


3 


1 





110) 


12 


3 


1 


1 





110) + |000) 


13 


5 


1 


1 





111) + 100) 


14 


7 


5 


1 





111) + 100) + |000) 


15 


11 


5 


1 






works in this case is the sp^ model.Bifl0 The basic fea- 
tures of this model are described in Table |, which lists 
the basis size (number of orbitals per atom) and the num- 
ber of independent Hamiltonian matrix elements for cou- 
pling between atoms in_the diamond structure out to 
third nearest neighbors O Table | also lists the proper- 
ties of other tight-bindmg models usedJji the literature. 



such as sp'^s* 



sp 



'rf^Mand 



sp 



^d'>s*m The number of 



parameters listed in this table is the number permitted 
by the symmetry of the model, which is not necessarily 
the same as that used in any specific implementation in 
the literature. 

For comparison, the bottom half of Table | lists the 
properties of several tight-binding models constructed 
from symmetrized coordinate eigenkets. The number of 
free parameters for the models generated by |100), |111), 
and 1 110) may be deduced easily from the link geome- 
try results presented in Appendix The corresponding 
numbers for the compound models (with more than one 
generator) were determined using the algorithm in Ap- 
pendix ^ 

The most striking feature of Table | is the relative 
paucity of free parameters in the symmetrized coordi- 
nate approach, which occurs because of the restriction 
(4.16) imposed by local gauge symmetry. Several of the 
symmetrized coordinate models are direct analogs of con- 
ventional models (i.e., they have identical symmetry); for 
example, the |111) (T^) model corresponds to sp^ , and 
t he 1 1 00) model corresponds to sp^d? [see Eqs. (2.10) and 
( ^.13| )]. However, the number of free parameters in con- 
ventional tight-binding theory grows steadily with dis- 
tance, whereas in the present theory there is no coupling 
beyond second nearest neighbors. 

This dearth of adjustable parameters means that the 
smallest basis sets do not provide a reliable model for 
the energy band structure. For example, in the [111) 
{Td) model, the splitting of the bonding and antibonding 
s states at the F point is the same as that of the p states 
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at r (and that of the p states at X). In the |100) model, 
there is no couphng at all between p orbitals on differ- 
ent atoms at the T point, so the splitting of bonding and 
antibonding states is zero. Such difficulties arise primar- 
ily because there is only one nearest-neighbor coupling 
parameter in these models. 

To increase the number of adjustable parameters with- 
out an undue increase in basis size, one must deliberately 
search for models with the most complicated topology 
available. A good starting point is the |111) {Oh) basis, 
which already has 3 nearest-neighbor parameters. Com- 
bining this with a |100) basis raises that number to 4 or 
5, depending on the relative values of the coordinates in 
the two generators. This 14-orbital model is a dramatic 
improvement over any of the smaller models; however, 
an extra |000) site adds substantial extra flexibility with- 
out much change in the basis size. Thus, the 15-orbital 
model generated by |000), |111), and |100) is probably the 
smallest basis capable of describing tetrahedral semicon- 
ductors accurately. In the language of conventional tight- 
binding theory, this would be referred to as an s^p^d^ f 
model. 

As shown in Table ||, this model has 17 free parameters 
(one of which is just the reference energy). Specific defi- 
nitions for these parameters are given in Table |l|, which 
also presents parameter values for Ge and Si obtained 
by fitting the band structure of the 15-orbital model to 
that given by the nonloc al . pinp irical pseudopotentials of 
Chelikowsky and Cohen.E3LjO Local Hamiltonian ma- 
trix elements are labeled V, whereas nonlocal on-site, 
nearest-neighbor, and second-nearest-neighbor terms are 
denoted a, f3, and 7, respectively. The subscripts a, 6, 
e, and / refer to independent sites generated by |0,0, 0), 
\r,r,r), \ — r, ~r, ~r) , and |r',0,0), respectively. The la- 
bels a, e, and / are the Wyckoff labels_£or sites of different 
symmetry in the diamond structure.Ej The label b is not 
correct Wyckoff notation; it actually refers to an inde- 
pendent e site, but the notation b was used here because 
these sites lie on the bonds between atoms. 

The energy band structure calculated from the param- 
eters in Table || is plotted in Fig. One can see that 
the 15-orbital model provides a good fit to the nonlocal 
pseudopotential bands from the bottom of the valence 
band to about 5 eV above the top of the valence band. 
Qualitative errors begin to occur near 9 eV at both the 
L and X points. For example, the Xi level near 9 eV in 
both figures should occur above 12 eV. This discrepancy 
can be eliminated, but the author has not found any way 
of doing so without adversely affecting the quality of the 
overall fit. 

It should be emphasized that the parameters in Ta- 
ble U are presented here merely as "proof of concept;" 
they are in no way intended as the final word on the sub- 
ject, and the author would be surprised if a better set 
were not found in the future. The quantities included in 
the fitting procedure were the valence- and conduction- 
band energy levels at F, X, L, and K. Effective masses 
and deformation potentials were not considered, and no 
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FIG. 2: Energy band structure of germanium and silicon. 
Solid lines: 15-orbital tight-binding model based on parame- 
ters in Table ^ Dotted lines: Nonlocal empirical pseudopo- 
tential model of Ref. The kinetic-energy cutoff for the 
latter calculation was (k -I- G)^ < 21(27r/ao)^, which corre- 
sponds to 113 plane waves at P. 



attempt was made at ensuring transferability. 

The main difficulty encountered during the fitting was 
the lack of any reliable method for establishing a sound 
starting point. Unlike the case for smaller tight-binding 
models, the present Hamiltonian has almost no simple 
analytical solutions (except for the F12', F12, and X2 
states, which are relatively unimportant) that can be 



used to determine starting values. The formula ( 4.42 ) for 
the kinetic energy provides a set of "free-particle" param- 
eters that is better than nothing, but in a 15-orbital basis, 
Eq. (4.42) is a poor approximation to the continuum ki- 
netic energy (5.19). After several months and dozens of 



different schemes (which still sampled only an infinitesi- 
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TABLE II: 


Independent Hamiltonian parameters in the 15-orbital model 


generated by |0,0, 0), \r,r,r), and 


lr',0,0). 




Parameter 




Value (Ry) 




Symbol 


Definition 


Ge 




Si 


Va 


{0,0,0|//|0,0,0) 


1.22536 




1.76282 




(r, r, r\H\r, r, r) 


1.18998 




1.65167 




{—r, r, r\H\ — r, r, r) 


1.07211 




1.17674 


Vf 


{r',0,OiiJjr',0,0) 


1.77902 




1.98485 


Ctab 


-{0,0,Oii/ir,r,r) 


-0.56483 




-0.50749 




-(0,0,Oi//i-r,r,r) 


0.07052 




0.24801 


aaf 


-(0,0,0|iJ|r-',0,0) 


0.06400 




-0.07980 




— (r, r, r\H\ — r, r, r) 


0.03841 




0.19402 


abf 


~lr,r,r\H\r' 
-~{-r,r,r\H\Q,r' ,Q) 


0.69832 




0.78464 




0.31958 




0.13684 




-\r',Q,Q\H\Q,r',Q) 


-0.43768 




-0.47311 


Pbb 


— {r, r, r\H\a — r, a — r, a — r) 


0.74306 




1.45766 


f3ee 


— {—r, r, r\H\a — r, a — r, a + r) 


-0.03011 




0.01537 


Pbf 


~{r,r,r\H\a ~ r' , a, a) 


-0.07041 




-0.35201 


Pef 


— {—r,r,r\H\a~r',a,a) 


0.28640 




0.40040 


Pff 


-{r',0,0\H\a,a-r',a) 


0.13140 




0.17401 


7ee 


— {— r, r, r\H\ — r,2a — r, 2a — r) 


0.02976 




0.02025 



mal fraction of parameter space), the author was unable 
to find any method whose success could honestly be at- 
tributed to anything other than trial and error. Hence, 
the development of a robust fitting procedure remains an 
unsolved problem. 



C. Dielectric function 

As a test of the field-particle coupling in the 15-orbital 
model, the imaginary part of the transiprse dielectric ten- 
sor was calculated from the formulacdO 



a/3 



(^) 



47r^e 



2^2 



E 



X 5{E, 



(27r)3 



(wk|p"|ck)(ck|/|t;k) 



(5.30) 



where hw is the photon energy, and |nk) is an eigen- 
ket of H with energy -Enk- The sum covered the four 
valence bands v and the seven lowest conduction bands 
c. The integral was perfornjed using a modified Gilat- 
Raubenheimer techniqucOoll3 based on 45961r-k points 
in the irreducible part of the Brillouin zondlsO (repre- 
senting 2048000 points in the full Brillouin zone). The 
energy interval for this calculation was 1 meV. 

To reveal more clearly the physical meaning of the cal- 
culated spectra, the same method was used to calculate 
the joint density of states function 



„ „, V ) J fix 



'fc. (5.31) 



The dielectric function differs from J{E) in that each 
transition is weighted by the oscillator strength 



Q/3 _ 



2{vk\p"\ck){ck\p'^\vk) 



(5.32) 



One may therefore use £2 and J to define an average 
oscillator strength at each energy by 



mwilo e2^(w) 



(5.33) 



In cubic crystals, the tensors (5.30) and (5.33) reduce to 



scalars: £3 (w) = t2{i^)5ap- 

The calculated dielectric function e2(cj) for Ge and Si 
IS plpi ted in Fig. |. This fi gure compares experimental 
datallS f or th e dielectric function with the values given 
by Eq. ( 5.30 ) for (i) the noxilQ^aLpseudopotential model 
of Chelikowsky and Cohen£jllJ'0 and (ii) the 15-orbital 
tight-binding model of Table ||. For each model, two 
plots of €2(0;) are given, corresponding to two different 
expressions for the momentum operator p. 

In pseudopotential calculations, optical properties a*a 
usually calculated from A • p coupling with p — — iW.O 
However, if the pseudopotential is nonlocal, this coupling 
is not gauge invariant; the correct linear coupling is give 
instead by the kinematic momentum p = (TO/i/i)[x, H]. 
Since pseudopotential calculations are usually performed 
in a plane- wave basis, a more convenient expression for 
the kinematic momentum is given by Eq. ( |5.9| ) (which is 
valid for both discrete and continuous coordinates x). 

In the present tight-bin ding theory, the kinematic mo- 
mentum is given by Eq. (5.2£). The two tight-binding 
functions plotted in Fig. H correspond to two different 
choices of the intra-atomic coordinates Xq , which are de- 
termined by the parameters r and r' in Table ||. One 
choice was the limit r ^ 0, r' ^ 0, which is equivalent 
to the zero-parameter model of Eqs. (1.3) and ( |1.4| ). The 
other, more physically realistic choice was r = and 
r' = ^/2r. The value r = was chosen because it yields 
equidistant lattice sites along the bond directions (111). 
The value r' — \/2r was used because a somewhat larger 
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FIG. 3: Imaginary part of the transverse dielectric func- 
tion of germanium and silicon. Dotted line: Experimen- 
tal data from Ref. |7|. Long dashed line: Nonlocal pseu- 
dopotential model of Ref. ^ with canonical (local) momen- 
tum p = —ihV. Dot-dashed line: Nonlocal pseudopoten- 
tial model of Ref. ^ with kinematic (nonlocal) momentum 
p = (m/ift)[x, //]. Short dashed line: 15-orbital tight-binding 
model from Table ^ with r — > and r' 0. Solid line: 15- 
orbital tight-binding model from Table ^ with r — and 
r' = V2r. 



value (e.g., 1.5r) breaks the link aff in Table |l|, whereas 
a somewhat smaller value (e.g., 1.3r) breaks the hnks 
f3bf and /?// (while simultaneously forming a new hnk 
Pbe)- These values of r and r' also generated successful 
starting values for some of the parameters in Table || 
(although the final fitted parameters were not very close 
to the starting values). 

Several conclusions may be drawn from Fig. ^. The 
first is that, within a given model (pseudopotential or 



tight binding), the choice of momentum operator does 
not have much numerical significance for the present cal- 
culation. This was to be expected on physical grounds, 
since the intra-atomic coordinate (in the tight-binding 
model) and the nonlocal part of the momentum (in the 
pseudopotential model) both lead to polarization effects 
within the atom. Yet it is well known that the bonds be- 
tween atoniSpare much easier to polarize than the atoms 
themselves J3'cl Hence, in a bulk semiconductor, intra- 
atomic effects yield only a minor numerical correction. 
This conclusion should remain valid in any system where 
the states are extended, but it may break-dpwn in sys- 
tems where localized states are important .l3E3 

The nonlocal part of the pseudopotential momentum 
tends to increase €2(1-0) in most frequency ranges, but it 
sometimes has the opposite effect (see, e.g., the region 
below 4 eV in Si). However, the intra-atomic coupling 
in the tight- binding model always decreases €2(1^-'). This 
may be understood by noting that the dominant nonlo- 
cal term in the Hamiltonian of Table || is the coupling 
f3bb along the bond between nearest neighbors. Increas- 
ing the value of r decreases the distance between b sites 
on neighborin g atom s, thereby decreasing the momentum 
matrix in Eq. ( 5.29| ) . This tends to increase (slightly) the 
discrepancy between the tight-binding and pseudopoten- 
tial dielectric functions. It is possible, however, that a 
different parametrization of the Hamiltonian might yield 
different results. 

The tight-binding and pseudopotential dielectric func- 
tions are quite similar in Ge, but there is a significant 
discrepancy at the E2 peak in Si. The reason for the 
difference between the models is apparent from the joint 
density of states J and average oscillator strength F plot- 
ted in Fig. ^. This figure shows that the tight-binding J 
is very accurate in Ge and somewhat less so in Si, as 
might have been expected from the quality of the fitted 
energy bands in Fig. |^. However, the tight-binding model 
underestimates the oscillator strength over almost the en- 
tire frequency range shown, typically by about 20%. The 
difference is most pronounced between 4 and 4.5 eV in 
Si, where it exceeds 30%. When combined with a slight 
underestimate of J in the same region, this leads to the 
discrepancy in the E2 peak noted above. 

The E2 peak in Si is associated with a volume in k 
space near (0.9, 0.1, 0.1)27r/ao,0 which is close to the X 
point. Thus, the physical reason for the error in the E2 
peak is probably the spurious Xi conduction band near 
9 eV in Fig. |[ Because this band is too low in energy, it 
mixes more strongly with the lowest Xi conduction band 
in the tight-binding model than it does in the pseudopo- 
tential model. This change in the wave function causes 
a corresponding change in oscillator strength. Thus, the 
majority of the error in the Si E2 peak would likely be 
eliminated if one could find an improved parameter set 
that raises the energy of the upper Xi conduction band. 

It is less clear whether the systematic underestimate 
of oscillator strength at all frequencies could be resolved 
by changing the Hamiltonian parameters. Oscillator 
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FIG. 4: Joint density of states J and average oscillator 
strength F for Ge and Si. Dotted line: Nonlocal pseudopo- 
tential model of Ref. |7^ with kinematic momentum p — 
{m/ih)]x, H]. Solid line: 15-orbital tight-binding model from 
Table y with r = |a and r' = ^/2r. 



strength was not included in the present fitting routine, 
so it is possible that with specific attention to this feature, 
one could improve the oscillator strength while maintain- 
ing the quality of the joint density of states. However, it 
is also possible that such an underestimate is a fundamen- 
tal limitation imposed by the small basis size in the tight- 
binding model. Thus, at present, the 15-orbital model is 
capable of providing semiquantitative predictions of os- 
cillator strength that reproduce all of the major trends 
exhibited by the pseudopotential model. Whether future 
developments bring it into precise quantitative agreement 
remains to be seen. 

Finally, it is worth noting that in Fig. ||, the calculated 
El peak (between 2 and 3 eV) for Ge is considerably less 
than the experimental value for both the pseudopoten- 



tial and tight-binding models. Chelikowsky and Cohen 
have attributed this discrepancy to the neglect of exci- 
ton effects £3 However, a recent first-principles calcula^. 
tion of 62 (ci^) that includes electron-hole interactionsEa 
shows that the contribution from excitons in Ge is not 
large enough to fill the gap, in F ig. ^. Thus, the empir- 
ical pseudopotential for GeESO probably needs further 
adjustment to increase the Ei peak. 



VI. DISCUSSION AND CONCLUSIONS 

This paper has shown that intra-atomic optical tran- 
sitions can be incorporated into tight-binding theory in 
a gauge-invariant way if the coordinate representation is 
taken as fundamental. Orthogonal atomic-like orbitals 
can be constructed from symmetrized coordinate eigen- 
kets, and the coupling to electromagnetic fields can then 
be described using lattice gauge theory. A model based 
on 15 such orbitals per atom is capable of describing the 
most important features of the tetrahedral semiconduc- 
tors Ge and Si. Tliifi basis is slightly larger than existing 
10-orbital modelsj^ but it has the advantages of (i) gauge 
invariance and (ii) providing an explicit wave function 
for the electron. A larger basis is needed in the present 
theory because the restrictions imposed by local gauge 
symmetry reduce the number of available Hamiltonian 
fitting parameters. 

The field-particle coupling derived bfTeiiji similar to 
that given by the Peierls substitution^BBEitaEa but 
the field-induced phase factor appears in the coordi- 
nate representation rather than the tight-binding repre- 
sentation. Thus, the present formalism includes intra- 
atomic coupling not present in|-.the Peierls substitution. 
Ismail-Beigi, Chang, and LouieB have recently presented 
a derivation of the Peierls phase for nonlocal Hamilto- 
nians in a continuous coordinate representation. They 
have argued that this derivation justifies the use of the 
Peierls substitution in tight-binding theory. However, 
their derivation cannot be extrapolated to tight-binding 
theory, because ordinary A • p coupling gives rise to intra- 
atomic interactions that are not included in the Peierls 
substitution. The existence of such interactions was not 
considered in the tight-binding theory of Ref. ||. 

It is interesting to consider whether there are any other 
ways of incorporating local gauge symmetry into tight- 
binding theory. One possibility is to work directly i n th e 
usual tight-binding representation [see, e.g., Eq. (1.4)], 



where the basis kets are labeled by the symmetry of the 
orbital and the position of the atom. If gauge symme- 
try is to be applied in this basis, the coordinate operator 
must be diagonal, hence all intr a-at omic matrix elements 
must be set to zero [as in Eq. (O)]. One could then in- 
troduce an Abelian U{1) gauge field on this lattice using 
the approach described above in Sec. [V, The results 
would be identical to those found in Sec. IV , except that 
the phase factor in the Hamiltonian (4.22ywould be ap- 
plied in the tight-binding representation rather than the 
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coordinate representation. Hence, this approach would 
constitute a "derivation" of the Peierls phase (|l.4|). Such 
a derivation would eli|iB|inate th e a mbiguity associated 
with the choice of pathcJ in Eq. (1.4). (Other techniques 
for eliminating path ambiguity are described in Refs. ^ 
and H.) 

The problem with this approach lies in its treatment of 
the coordinate operator. In the theory described in Sec. 
IV , when the basis size is increased, the eigenvalue spec- 
trum of the coordinate operator remains nondegenerate, 
tending (in the limit of infinite basis dimensions) toward 
a continuous spectrum. However, if the coordinate oper- 
ator is required to be diagonal in the tight-binding ba- 
sis, its eigenvalue spectrum is always degenerate, tend- 
ing (in the limit of infinite basis dimensions) toward a 
discrete spectrum with infinite degeneracy. Hence, any 
tight-binding theory that is either based on or equivalent 
to the Peierls substitution cannot reproduce the correct 
continuum limit of the coordinate operator. 

As a generalization of the above approach, one 
migbLp^Jstip-cpMpj^^ a non-Abelian gauge 

fieldBeemyy-BBEiy in the tight-binding basis. 
The idea would be to treat the tight-binding electron 
as a new type of "elementary particle" with some in- 
ternal degrees of freedom (corresponding to the symme- 
try labels of the atomic orbitals) that are coupled to the 
gauge field. In this way, one might hope to reproduce 
the effects of intra-atomic coupling while remaining in 
the tight-binding basis. There are, however, numerous 
difficulties with this approach. 

First, the lattice sites in lattice gauge theory represent 
states of the same particle at different positions. Hence, 
these states are identical apart from their positions. How- 
ever, in tight-binding theory the atoms are generally not 
the same. Second, the field equations for a non-Abelian 
gauge field are intrinsically nonlinear, because the field 
carries its own charge and is coupled directly to itself (i.e., 
it is self-radiating) . It is therefore difficult to imagine how 
such a field could reproduce ordinary electromagnetism, 
in which the field has no charge, and nonlinearities arise 
only from interactions with matter. Third, one would 
need to define a new gauge field theory every time one 
added new orbitals to the model, and every one of these 
non-Abelian field theories would need to reproduce the 
results of Abelian electromagnetic theory. Finally, the 
coordinate operator in this approach would still have a 
discrete, degenerate eigenvalue spectrum. 

Thus, it appears that the present approach — that 
is, an Abelian U{1) gauge field in the coordinate 
representation — is the only gauge-invariant method for 
including electromagnetic fields in empirical tight- 
binding theory that tends toward the correct continuum 
limit as the basis dimensions are increased. In this case, 
the only way that the essential structure of the theory 
can be modified is to change the topology of the sys- 
tem so as to increase the number of links between lattice 
sites. This would increase the number of free parameters 
in the Hamiltonian, thereby permitting a reduction in 



basis size. Such a modification would clearly be benefi- 
cial, but it is not obvious that there exists any alternative 
topology for general lattices that is capable of reproduc- 
ing continuum electromagnetism unambiguously. Hence, 
this possibility will not be explored further here. 
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APPENDIX A: SYMMETRIZED ORBITALS 

This section presents symmetrized orbitals obtained by 
applying the symmetry operations of the cubic group Oh 
to the coordinate eigenkets |100) and |110). [The orbitals 
for |111> may be found in Eq. (|2.1l|) .] For |100), if_the 
basis kets are ordered as {|100), |010), |001), |100), |010), 
1 001)}, then the symmetrized orbitals are 



iri) 



|r?2) 



75(1,1,1,1,1,1), 

^(0,0,1,0,0,-1), 
^(-1,-1,2,-1,-1,2), 

i(i, -1,0,1, -1,0). 



(Al) 



For |110), if the basis kets are ordered as {|011), |011), 
1011), 1011), |101), |101), |101), |101), |110), |110), |110), 
1 110)}, then the symmetrized orbitals are 



iri) = 
ir^s) - 
K2) = 



|r25'': 



25/ 



= i7^(l, 1,1, 1,1, 1,1, 1,1, 1,1,1), 

= ^(1,-1,1,-1,1,-1,1,-1,0,0,0,0), 

= ^(-1,-1,-1,-1,-1,-1,-1,-1,2,2,2,2), 

= -1^(1,1,1,1,-1,-1,-1,-1,0,0,0,0), 

= i(0,0,0,0,0,0,0,0,l, -1,-1,1), 

= l/z(x2-y2)) 

= i7i(-l,l,-l,l,l,-l,l,-l,0,0,0,0). (A2) 



For the triply degenerate representations Fis, r25', and 
r25, only one representative orbital is given; the others 
may be obtained from cyclic permutations of a:, y, and z. 
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APPENDIX B: GEOMETRY OF VORONOI 
POLYHEDRA 

This appendix presents an algorithm for calculating 
the geometry of the Voronoi polyhedra associated with a 
given set of nodes . The basic element in this algorithm 
is a procedure for finding the edges of the polyhedra. 
An edge of a Voronoi polyhedron is a finite line segment 
consisting of points that are closer to three (or more) 
nodes than to any other nodes. The first step is therefore 
to determine the equation defining this line. 

Any three noncoUinear points Xj, Xj, and x^, define a 
plane whose normal is the vector 



X7 X X 7 



-X, X Xfc +Xfe X Xi, 



(Bl) 



where dji = xj — x^. This plane is the set of points x 
satisfying 



riy fe • (x - Xi 



0. 



(B2) 



The line consisting of all points x equidistant from Xj, 
Xj, and Xfc may therefore be written as 



X = Xyfe + Xh, 



(B3) 



where A is a real parameter, xiijk = riyfe /riyfc, and x^fe 
is the point in the plane (B2) equidistant from x^, Xj, 
and Xfc. To determine this point, note that points x 
equidistant from x^ and Xj satisfy 



(Xj - X.,) • [X - i(Xj + Xi)] = 0. 



(B4) 



The point x^j/c therefore satisfies the three equations 



(r = l,2,3) 



(B5) 



in which 



ai = d,i. 



a2 



^3 — T^ijk, 



(B6) 



which may be viewed as a set of oblique (more specifi- 
cally, monoclinic) basis vectors, and 



1 



•^1 

= 'x.i ■ (xj X Xfc). 



(B7) 



The solution to Eqs. (|B5|) is given by 



(B8) 



s=l 



where bs is a reciprocal basis vector satisfying a,. • b<, = 
(5r.s; e.g.. 



bl 



a2 X as 



(B9) 



ai • (a2 X as) 

Now since ai • (a2 x a^) — nfji,, the vectors b^ are given 
explicitly by 



bl = [dji{dl,) - dki{dj, ■ dk^)]/^ 
ba = [dk^{(fji) - dji{dj, ■ dkt)]/r 
bg = n,jk/njjk. 



2 

'ijk 1 



(BIO) 



metric form 



Equation (B8) may then be rearranged in the more sym- 



^ijk 



XiM?fe(dji • (ik^)] + Xj[dl^{dkj ■ dij)] + Xfc[d?-(djfc • djk)] 



(Bll) 



r 



This result, together with Eq. (B3), defines the line 
equidistant from nodes x^ , Xj , and x^ . 

The next step is to determine whether any segment of 
this line forms an edge of a Voronoi polyhedron. Points 
on such a segment must lie closer to x^, Xj, and x^ than 
to any other node x;. For each node xj, one calculates 



il ' '^ijk ■ 



(B12) 



If ai — 0, then x; lies in the plane (B2). In this case, if 
|xi — Xijk\ < |xi — Xyfcl, then no portion of the line (B3) 



forms an edge of a Voronoi polyhedron. On the other 
hand, if jx; — Xy^l > jx^ — Xy^l, then the line (B3) may 



form an edge (depending on the position of the other 
nodes xj'). 



If ai 7^ 0, then x; does not lie in the plane (B2). In 
this case, points on the line ( p^ ) that are closer to x^ 
than to x; satisfy [cf. Eq. (B4)] 



dii ■ (xyfe + Xhijk - Xij) > 0, (B13) 

in which x^; — ^{xi + x;). Hence, the position of the 
point equidistant from x^ , Xj , x^ , and x; i s giv en by the 
following value of the parameter A in Eq. (B3): 

dii ■ (xiz — Xijk) 



A; 



(B14) 
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FIG. 5: Partitioning of the surface Sij 
area calculation in Eq. (B2C). 



into triangles for the 



One may then define 

Amin = max(A/|a; > 0) 



(B15) 



(i.e., the maximum value of A; for all I such that ai > 0) 
and 



An 



min(Ai|ai < 0). 



(B16) 



If Aniax > Amin, then the line segment (B3) with Amin < 
A < Amax forms an edge of a Voronoi polyhedron. This 
establishes the positions of two corners of the polyhedron: 



Xc 
Xf,' 



^ijk Amin Aijfc 5 



(B17) 



The set of all nodes in the plane (B2) that lie closer 
to the line segment Amin < A < Amax than any other 
node defines what is called a plaquette. Since there are 
in general more than three such nodes, it is convenient 
to define a unique label q for each plaquette, with 



(B18) 



for any members x^, Xj, and x^; of the given plaquette. 
(The sign of is fixed by some convention for the or- 
dering of the nodes x^, x^, and x^.) Each plaquette is 
associated uniquely with one edge of a Voronoi polyhe- 
dron, the length of which is 

dq = |Xc' — Xc| = Amax ^ Amin, (-^19) 

with > by definition. 

At this point, one has sufficient information to deter- 
mine whether a link exists between any pair of nodes x.^ 
and Xj. The first step is to use the above procedure to 
find all of the corner points Xc common to nodes i an d 
j. By definition, all such points lie in the plane (B4). 



The set of these points defines a polygon, the perimeter 
of which consists of the line segments ( B19 ). The area of 
the polygon may be calculated by numbering the corner 
points Xc in sequential order around the perimeter of the 
polygon, then partitioning the polygon into triangles as 
shown in Fig. 0. The normal area vector is 



Xi) X (Xc - Xi), 



(B20) 



where Nij is the number of corner points common to 



nodes 



is the area of the 



I and j. The area Sij — f^^-^ 
surface shared by the Voronoi polyhedra for sites i and 
j; note that if Nij — 1 or 2, the shared region is a point 
or line, and the area (B2C) is zero. Nodes x^ and Xj are 
linked only if Sij > 0. 

The volume fli of the Voronoi polyhedron for node Xj 
may be calculated from Sij and dij. One simply inte- 
grates the identity V • x = 3 over the polyhedron, using 
the divergence theorem and the fact that the plane con- 
taining Sij is the perpendicular bisector of dij (although 
dij need not intersect Sij itself). The result is 



(B21) 



For each link £ = one can construct a polyhedron 

by drawing lines from nodes x^ and x^ to each of their 
common corner points Xc. The volume of this polyhedron 
is 



1 



:Spdi. 



(B22) 



The volume Vli = flij is bisected by Se = Sij, with half 



lying in Q,i and half in Vlj ; hence 



(B23) 



The nodes in plaquette q define a polygon in the plane 



( B2 ) ; the perimeter of this polygon is formed by the links 
dg. Hence, the area of the plaq uette can be calculated in 
the same way as the link area ( B20| ): 



-Y 

2 ^ 

4=3 



(xi-i - xi) X (xi - xi). 



(B24) 



where Nq is the number of nodes in plaquette q. A poly- 
hedron may be constructed for each plaquette by drawing 
l ines from each node of the plaquette to the corner points 
(B17); the volume of this polyhedron is 



^Sqdq. 



(B25) 



Finally, the plaquette surfaces Sq partition all of space 
into nonoverlapping polyhedra (this is referred to as a 
Delaunay tessellatioiiE3). These polyhedra (or cells) are 
in one-to-one correspondence with the corner points Xc 
of the Voronoi polyhedra. The volume of cell c is 



(B26) 



q€c 



where the direction of Sq is chosen to point outward |fcom 
flc (note that Xc does not necessarily lie inside flcf^ so 
the dot product may be negative for some q). 
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A useful set of sum rules for verifying the consistency 
of a calculated geometry is 



(B27) 



where is the volume of some region over which the 
node distribution is periodic, such as a primitive cell in 
a Bravais lattice (not to be confused with the generally 
nonperiodic cell Oc)- The sum rule for Qi follows directly 
from the definition of fli given in Sec. [II, since every 



point in f2 must lie in at least one Voronoi polyhedron, 
and the only regions of overlap between polyhedra are 
points, lines, or planes of zero volume. The sum rule for 
ric was proven in Ref. 1^. The sum rule for fli follows 
from that for 17^, since the set {fif} is jus t another way of 
partitioning the set {fli} [see Eq. ( B23| )]. Likewise, the 



sum rule for flq follows from that for f2c, since the set 
{Oq} is jus t another way of p artitioning the set {ftc} [see 
Eqs. pl|), (HI), and (p26|)]. 



APPENDIX C: EXAMPLES OF LINK 
GEOMETRY 

This appendix presents values of the link lengths di 
and surface areas for several lattices. The simplest 
geometry occurs for Bravais lattices, of which only the 
cubic lattices are considered here. For the simple cubic 
lattice, only nearest neighbors are linked, with di = a 



and S'l 



where a is the lattice constant. For the 



body-centered cubic lattice, both first and second nearest 
neighbors are linked, with 



and 



V3 3\/3 2 

z lb 



d2 = a, 82 = -a^ 



(CI) 



(C2) 



For the face-centered cubic lattice, only nearest neighbors 
are linked, with 



di = 



^/2' 



Si = 



4V2' 



(C3) 



The remaining lattices to be considered are those ob- 
tained by putting symmetrized orbitals on the atomic 
sites of the diamond or zinc-blende structure. If only a 
single s orbital per atom is used (i.e., one |000) basis ket 
per atom), then each atom is linked to 4 nearest neigh- 
bors and 12 second-nearest neighbors, with 



and 



di = \/3a, 5*1 = sVSa^, 
/2 

d2 = 2V2a, S2 = 



(C4) 



(C5) 



Here a = jUq, where qq is the conventional cubic lattice 
constant. Note that in this case, the link ^2 does not 
intersect the surface 82- 

Coupling between second-nearest neighbors persists in 
models with more than one orbital per atom. In the 
"sp^" model with four [111) sites per atom (generated 
by applying the symmetry operations of to \r,r,r) 
and \a — a — r, a — r)) , each site is linked to three others 
on the same atom: 

do = 2V2r, So = (C6) 



one on a neighboring atom: 

di = \/3(a - 2r), Si^sVSa^, (C7) 
and three second-nearest neighbors: 



d2 = 2V2{a-r), 82 = ^^?. 



(C8) 



In the model generated by either Td or Oh and |r, 0, 0), 
each of the six sites is linked to four others on the same 
atom: 



do = \/2r, 5*0 = 



a^(5a — 6r) 
4V2(a - r) ' 



four nearest neighbors: 



2 7 

di = V3a2-4ar + 2r2, " ^ 



2(a - r) 



and four second-nearest neighbors: 



d2 = V2(2a-r), ^2 = 



a^{a - 2r) 
4V2(a - r) 



(C9) 



(CIO) 



(Cll) 



In the model generated by either or Oh and |r, r, 0), 
each of the 12 sites is linked to four others on the same 
atom, two of which have 

do = V2r, So = J^^^, (C12) 

and two of which have 

rfo = ^/2., = (C13) 

4v2(a — r) 

Each site is also linked to two nearest neighbors: 

a^di 



di = 3a^ — 8ar + 6r^, Si 



and one second-nearest neighbor: 



d2 = 2\/2(a-r), ^2 



2{a~r)' 



2V2{a - r) ' 



(C14) 



(C15) 
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In the model generated by Oh and \r,r,r), there are 
two distinct lattice sites. Sites such as \r,r,r) and \a — 
r,a — r,a — r) are labeled b because they lie on the bonds 
between atoms, whereas sites such as | — r, — r, — r) and 
\a+r, a+r, a+r) are labeled e because they lie on "empty" 
bonds. (Both of these sites are actually WyckofF e sites, 
but they are inequivalent, because the site symmetry of 
the atoms in diamond is Tj,.) Each b site is linked to 
three e sites on the same atom (and vice versa), with 



nbe 



a{9a - 4r) 



(C16) 



Each b site is linked to one nearest-neighbor b site: 

df = V3(a - 2r), 5f = (C17) 



and three nearest-neighbor e sites: 



< = v/a2 + 2(a-2r)2, fi-jx^ = -a<. (C18) 



Each e site is also linked to three b sites via (CIS), to six 
nearest-neighbor e sites via 



1 



= v/2a2 + (a-2r)2, = ^adf, (C19) 

and to three second-nearest neighbor e sites via 

/2 

df = 2V2(a-r), 5*^ = — a(a -|- 2r). (C20) 
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